/**
 * file name  : Jenoware/it/ciano/jenoware/Haplotype.java
 * authors    : Luciano Xumerle
 * created    : mar 23 ott 2007 16:37:16 CEST
 *
 * Copyright (c) 2004-2005 Luciano Xumerle. All rights reserved.
 *
* This file is part of gephinfo.
*
* gephinfo is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program.  If not, see <http://www.gnu.org/licenses/>.

 *
 */

package it.ciano.jenoware;

import java.util.ArrayList;
import java.util.Collections;
import java.util.Arrays;
import it.ciano.util.SpecialFunction;

/**
 * The class extends Value and stores the allele value (String or numeric).
 *
 * @author Luciano Xumerle
 * @version 0.0.1
 */
final public class Haplotype
{

    // the haplotype is an array of Genotype objects.
    // Each haplotype filament is defined by the alleles' order of each Genotype.
    private Genotype[] _haplo;

    /**
     * Class Constructor.
     */
    public Haplotype ()
    {
        _haplo=null;
    }

    /**
     * Class Constructor.
     *
     * @param haplotype array of map ordered genotypes.
     */
    public Haplotype ( Genotype[] haplotype )
    {
        _haplo=haplotype;
    }

    /**
     * Class Constructor.
     *
     * @param haplotype the PairAllele list.
     */
    public Haplotype ( PairAllele[] haplotype )
    {
        _haplo=new Genotype[haplotype.length];
        for ( int i=0; i<haplotype.length; i++ )
        {
            String a="Marker"+i;
            _haplo[i]=new Genotype ( a, haplotype[i].getAllele ( 1 ),
                    haplotype[i].getAllele ( 2 ), "", 1.0,
                    haplotype[i].isMantainAllelesOrder() );
        }
    }

    /**
     * Returns the number of markers of the haplotype.
     *
     * @return the number of markers of the haplotype.
     */
    public int size()
    {
        if ( _haplo == null )
            return 0;
        return _haplo.length;
    }

    /**
     * Returns the Genotype object in given position.
     * The position values are from 0 to getMarkersNumber()-1
     *
     * @param i Genotype position inside haplotype.
     * @return the Genotype object in given position.
     */
    public Genotype getGenotypeInPosition ( int i )
    {
        if ( i>-1 && i<size() )
            return _haplo[i];
        return null;
    }

    /**
     * Returns the PairAllele object composed by the haplotype.
     * The PairAllele alleles' strings are simply the concatenation between alleles
     * in each markers.
     *
     * @return the PairAllele object.
     */
    public PairAllele toPairAllele()
    {
        StringBuffer all1=new StringBuffer();
        StringBuffer all2=new StringBuffer();
        for ( int i=0; i<_haplo.length; i++ )
        {
            if ( _haplo[i].isValid() )
            {
                all1.append ( _haplo[i].getAll1() );
                all2.append ( _haplo[i].getAll2() );
            }
            else
            {
                all1.append ( '-' );
                all2.append ( '-' );
            }
        }
        return new PairAllele ( all1.toString(), all2.toString() );
    }


    /**
     * Returns the Haplotype formatted as Gerbil input file format.
     * Each line represents a haplotype.
     * Each two consequative lines correspond to one person.
     * Each SNP is 0,1,2 or ? separated by a tab.
     * A heterozygous site that its phasing is unkown is represented by '2'.
     * A missing data entry is represented by '?'.
     *
     * @param firstAllele The list of first alleles for each SNP in Haplotype.
     * @return The array containing the two consequative lines corresponding to one person.
     */
    public String[] toGerbil ( String[][] firstAllele )
    {
        if ( _haplo.length != firstAllele.length )
            return null;
        StringBuffer all1=new StringBuffer();
        StringBuffer all2=new StringBuffer();
        for ( int i=0; i<_haplo.length; i++ )
        {
            if ( _haplo[i].isValid() )
            {
                if ( !_haplo[i].isHomozigote() )
                {
                    all1.append ( '2' );
                    all2.append ( '2' );
                }
                else if ( _haplo[i].getAll1().equals ( firstAllele[i][0] ) )
                {
                    all1.append ( '0' );
                    all2.append ( '0' );
                }
                else
                {
                    all1.append ( '1' );
                    all2.append ( '1' );
                }
            }
            else
            {
                all1.append ( '?' );
                all2.append ( '?' );
            }
            if ( i+1<_haplo.length )
            {
                all1.append ( "\t" );
                all2.append ( "\t" );
            }
        }
        return new String[] {all1.toString(), all2.toString() };
    }


    /**
     * Returns the Haplotype formatted as Gerbil input file format.
     *
     * @param firstAllele The list of first alleles for each SNP in Haplotype.
     * @return The String containing the line corresponding to one person haplotype.
     */
    public String toSNPhap ( String[][] firstAllele )
    {
        if ( _haplo.length != firstAllele.length )
            return null;
        StringBuffer all=new StringBuffer();
        for ( int i=0; i<_haplo.length; i++ )
        {
            if ( _haplo[i].isValid() )
            {
                if ( !_haplo[i].isHomozigote() )
                {
                    all.append ( '1' );
                    all.append ( "\t" );
                    all.append ( '2' );
                }
                else if ( _haplo[i].getAll1().equals ( firstAllele[i][0] ) )
                {
                    all.append ( '1' );
                    all.append ( "\t" );
                    all.append ( '1' );
                }
                else
                {
                    all.append ( '2' );
                    all.append ( "\t" );
                    all.append ( '2' );
                }
            }
            else
            {
                all.append ( '0' );
                all.append ( "\t" );
                all.append ( '0' );
            }
            if ( i+1<_haplo.length )
            {
                all.append ( "\t" );
            }
        }
        return all.toString();
    }


    /**
     * Returns the Haplotype formatted as HIT input file format.
     * ? means missing value.
     *
     * @param firstAllele The list of first alleles for each SNP in Haplotype.
     * @return The array containing the two consequative lines corresponding to one person.
     */
    public String[] toHIT ( String[][] firstAllele )
    {
        if ( _haplo.length != firstAllele.length )
            return null;
        StringBuffer all1=new StringBuffer();
        StringBuffer all2=new StringBuffer();
        for ( int i=0; i<_haplo.length; i++ )
        {
            if ( _haplo[i].isValid() )
            {

                if ( _haplo[i].getAll1().equals ( firstAllele[i][0] ) )
                    all1.append ( '1' );
                else
                    all1.append ( '2' );

                if ( _haplo[i].getAll2().equals ( firstAllele[i][0] ) )
                    all2.append ( '1' );
                else
                    all2.append ( '2' );
            }
            else
            {
                all1.append ( '?' );
                all2.append ( '?' );
            }
        }
        return new String[] {all1.toString(), all2.toString() };
    }

    /**
     * Returns the Haplotype formatted as HAP input file format.
     * Each line represent an individual (id &lt;TAB&gt; haplotype).
     * 'A','G','C' or 'T' representing a homozygous genotype.
     * 'H' representing a heterozygous genotype.
     * '?' representing a missing genotype.
     *
     * @param destLabel The array containing the first allele label and the labels [A,C,G,T] to be used.
     * @return the line representing an individual.
     */
    public String toHAP ( String[][] destLabel )
    {
        StringBuffer all=new StringBuffer();
        for ( int i=0; i<_haplo.length; i++ )
        {
            if ( _haplo[i].isValid() )
                if ( _haplo[i].isHomozigote() )
                    if ( _haplo[i].getAll1().equals ( destLabel[i][0] ) )
                        all.append ( destLabel[i][1] );
                    else
                        all.append ( destLabel[i][2] );
                else
                    all.append ( 'H' );
            else
                all.append ( '?' );
        }
        return all.toString();
    }

    /**
     * Returns the list of PairAllele objects containing every possible haplotype pairs
     * observable using the given Genotype list.
     *
     * @param in the ordered list of Genotype objects.
     * @return the list of possible PairAllele. Each allele of each PairAllele is a haplotype.
     */
    public static PairAllele[] snp2haplo ( Genotype[] in )
    {
        ArrayList res=new ArrayList();
        int[][] phases=SpecialFunction.getCombinationList ( in.length );

        for ( int j=0;j<phases.length;j++ )
        {
            // combinations starting with 1 are mirrored to
            // combinations starting with 0.
            if ( phases[j][0]==1 )
                continue;
            Genotype[] aa=new Genotype[in.length];
            for ( int i=0;i<in.length;i++ )
            {
                aa[i]=new Genotype ( in[i].getName(),in[i].getAllele ( 1 ),
                        in[i].getAllele ( 2 ), "", 1.0,
                        ( phases[j][i]==1 ) );
            }
            res.add ( new Haplotype ( aa ).toPairAllele() );
        }
        return getUniquePairAlleleValue ( res );
    }


    /**
     * Returns the list of PairAllele objects containing every possible haplotype pairs
     * observable using the given Genotype list.
     *
     * @param in the ordered list of Genotype objects.
     * @param haplotypes The available haplotypes.
     * @return the list of possible PairAllele. Each allele of each PairAllele is a haplotype.
     */
    public static PairAllele[] snp2haplo ( Genotype[] in, MarkerFreq haplotypes )
    {
        Allele[] list=haplotypes.getOrderedAlleleList();

        ArrayList res=new ArrayList();

        for ( int j=0; j<list.length; j++ )
        {
            StringBuffer second=new StringBuffer();
            String haplo=list[j].toString();
            boolean ok=true;
            for ( int i=0; i<in.length; i++ )
            {
                final String haplochar=haplo.substring ( i, i+1 );
                if ( in[i].getAll1().equals ( haplochar ) )
                    second.append ( in[i].getAll2() );
                else if ( in[i].getAll2().equals ( haplochar ) )
                    second.append ( in[i].getAll1() );
                else
                {
                    ok=false;
                    break;
                }

            }
            if ( ok )
            {
                PairAllele tt=new PairAllele ( haplo, second.toString() );
                if ( tt !=null )
                {
                    tt.setFrequency ( haplotypes.getUnrelatedGenotypeFreq ( tt.getConcatAlleles() ) );
                    res.add ( tt );
                }
            }
        }
        return getUniquePairAlleleValue ( res );
    }

    final private static PairAllele[] getUniquePairAlleleValue ( ArrayList res  )
    {
        ArrayList rr=new ArrayList();
        for ( int i=0;i<res.size();i++ )
        {
            PairAllele now= ( PairAllele ) res.get ( i );
            boolean missing=true;
            for ( int j=0;j<rr.size();j++ )
                if ( now.equals ( ( PairAllele ) rr.get ( j ) ) )
                    missing=false;
            if ( missing )
                rr.add ( now );
        }
        return ( PairAllele[] ) rr.toArray ( new PairAllele[rr.size() ] );
    }

    /**
     * Overriding toString method.
     *
     * @return the Allele value.
     */
    public String toString()
    {
        StringBuffer a=new StringBuffer();
        for ( int i=0; i<_haplo.length; i++ )
        {
            a.append ( _haplo[i].getConcatAlleles() );
            a.append ( "\n" );
        }
        return a.toString();
    }

    /**
     * Overriding equals method.
     *
     * @param o the object to compare.
     * @return true or false.
     */
    public boolean equals ( Object o )
    {
        Haplotype a= ( Haplotype ) o;
        for ( int i=0; i<_haplo.length; i++ )
            if ( !_haplo[i].equals ( a.getGenotypeInPosition ( i ) ) )
                return false;
        return true;
    }

    /**
     * Overriding clone method.
     *
     * @return a new cloned object ID
     */
    public Object clone()
    {
        Genotype[] a=new Genotype[_haplo.length];
        for ( int i=0; i<_haplo.length; i++ )
            a[i]= ( Genotype ) _haplo[i].clone();
        Haplotype cloned = new Haplotype ( a );
        return ( Object ) cloned;
    }

    /**
     * Returns the value of D.
     * A deviation of the observed frequencies from the expected
     * is referred to as the linkage disequilibrium parameter,
     * and is commonly denoted by a capital D [3] as defined by:
     * <i>D = x11 - p1q1</i>.
     * <i>x11</i> is the observed frequency of the haplotype
     * and <i>p1q1</i> is the expected frequency.
     *
     * @param a1 The first marker observed frequencies of alleles.
     * @param b1 The second marker observed frequencies of alleles.
     * @param haploFreq The observed haplotype frequency.
     * @return the value of D.
     */
    public static double getD ( double a1, double b1, double haploFreq )
    {
        return  ( haploFreq - ( a1*b1 ) );
    }

    /**
     * Returns the value of D'.
     * <i>D'=D/D<sub>max</sub></i> if <i>D &gt;= 0</i>.
     * <i>D'=D/D<sub>min</sub></i> if <i>D &lt; 0</i>.
     * D<sub>max</sub> is given by the smaller of <i>a1b2</i> and <i>a2b1</i>.
     * D<sub>min</sub> is given by the larger of <i>-a1b1</i> and <i>-a2b2</i>.
     *
     * @param D The D value.
     * @param a1 the frequency of the first allele in the first marker.
     * @param a2 the frequency of the second allele in the first marker.
     * @param b1 the frequency of the first allele in the second marker.
     * @param b2 the frequency of the second allele in the second marker.
     * @return Returns the value of D'.
     */
    public static double getD1 ( double D, double a1, double a2, double b1, double b2 )
    {
        if ( D>=0 )
        {
            return ( D / ( Math.min ( a1*b2, a2*b1 ) ) );
        }
        else
        {
            return D / ( Math.max ( -a1*b1, -a2*b2 ) );
        }
    }


    /**
     * Returns the value of R<sup>2</sup>.
     * Denoted as <i>R<sup>2</sup>= D<sup>2</sup> / a1*a2*b1*b2</i>.
     *
     * @return the value of R<sup>2</sup>.
     */
    public static double getR2 ( double D, double a1, double a2, double b1, double b2 )
    {
        return ( D*D ) / ( a1*a2*b1*b2 );
    }


    /**
     * Returns an error rate detected between two strings.
     * It is calculated as (number of unmatched chars)/(total number of char).
     *
     * @param expected the expected string value.
     * @param observed the observed string value.
     * @return the error rate.
     */
    public static double getSimilarity ( String expected, String observed )
    {
        char[] ex = expected.toCharArray();
        char[] ob = observed.toCharArray();
        int min=Math.min ( expected.length(), observed.length() );
        int max=Math.max ( expected.length(), observed.length() );
        double error=max-min;
        for ( int i=0;i<min;i++ )
            if ( ex[i]!=ob[i] ) error+=1;
        return error/ ( double ) max;
    }


    /**
     * Returns the switch accuracy value calculated with the Shin Lin method.
     * Article: Haplotype inference in random population samples (2002).
     *
     * @param ex1 the first expected haplotype.
     * @param ex2 the second expected haplotype.
     * @param ob1 the first observed haplotype.
     * @param ob2 the first observed haplotype.
     * @return the switch accuracy value.
     */
    public static double getSwitchAccuracy ( String ex1, String ex2, String ob1, String ob2 )
    {
        String[] ex=new String[] {ex1, ex2};
        String[] ob=new String[] {ob1, ob2};
        Arrays.sort ( ex );
        Arrays.sort ( ob );
        System.out.println ( "or1: " + ex[0] );
        System.out.println ( "or2: " + ex[1] );
        System.out.println ( "in1: " + ob[0] );
        System.out.println ( "in2: " + ob[1] );
        char[] e1 = ex[0].toCharArray();
        char[] e2 = ex[1].toCharArray();
        char[] o1 = ob[0].toCharArray();
        char[] o2 = ob[1].toCharArray();
        int countSwitch=0;
        int countHET=0;
        boolean isSwitched=false;
        for ( int i=0; i<e1.length; i++ )
        {
            if ( e1[i]==e2[i]  )
            {
                if ( e1[i]!=o1[i] || e1[i]!=o2[i] )
                    System.err.println ( "ERROR IN HAPLOTYPE (getSwitchAccuracy)" );
                continue;
            }
            countHET++;
            if (  !isSwitched && e1[i]==o2[i] && e2[i]==o1[i] )
            {
                isSwitched=true;
                countSwitch++;
            }
            if (  isSwitched && e1[i]==o1[i] && e1[i]==o1[i] )
            {
                isSwitched=false;
                countSwitch++;
            }
        }
        System.out.println ( "countSwitch: " + countSwitch + " countHET: " + countHET );
        if ( countSwitch==0 )
            return 1;
        return  ( double ) ( countHET -1-countSwitch ) / ( double ) ( countHET -1 );
    }


    /**
     * Returns the list of compatible haplotypes stored as PairAllele objects.
     * Each PairAllele object contains two haplotype and the probability of the
     * subject to have the two haplotypes. The first PairAllele object has the
     * higher probability.
     *
     * @param sbj the subject.
     * @param snpsList list of SNP's names to be used.
     * @param haploFreq the frequency of each haplotipe.
     * @return the list of compatible haplotypes stored as PairAllele objects.
     */
    public static PairAllele[] inferPossibleHaplotypes ( Sbj sbj, String[] snpsList, MarkerFreq haploFreq )
    {
        /**
        * gets the genotype list used to get haplotypes
        */
        Genotype[] ain = new Genotype[snpsList.length];
        for ( int i=0; i<snpsList.length; i++ )
            ain[i]=sbj.getGenotype ( snpsList[i] );

        /**
        * The list of Genotypes is used to get the list
        * of each possible haplotype pairs.
        */
        PairAllele[] ahaplo= Haplotype.snp2haplo ( ain, haploFreq );

        /**
        * for each haplotype pair the program gets the frequency of
        * the two haplotypes using the haploFreq object.
        * The sum is the sum of each frequency used.
        * The sum is used to get the probability (frequency/sum) to observe the two
        * haplotype in the subject.
        * The second line of array contains the frequency/sum values.
        */
        double sum=0;
        double[] afreqs=new double[ahaplo.length];
        for ( int aa=0; aa<ahaplo.length; aa++ )
        {
            afreqs[aa]=haploFreq.getUnrelatedGenotypeFreq ( ahaplo[aa].getConcatAlleles() );
            sum+=afreqs[aa];
        }
        for ( int aa=0; aa<ahaplo.length; aa++ )
        {
            afreqs[aa]=afreqs[aa]/sum;
            ahaplo[aa].setFrequency ( afreqs[aa] );
        }

        /**
        * sort the frequencies and reorder the PairAllele list to
        * get and ordered list where the first element has
        * the higher probability.
        */
        Arrays.sort ( afreqs );
        PairAllele[] res=new PairAllele[ahaplo.length];
        int resind=afreqs.length-1;
        for ( int i=0; i<afreqs.length; i++ )
            for ( int a=0; a<ahaplo.length; a++ )
            {
                if ( ahaplo[a]==null ) continue;
                if ( ahaplo[a].getFrequency() == afreqs[i] )
                {
                    res[resind]=ahaplo[a];
                    resind--;
                    ahaplo[a]=null;
                }
            }
        return res;
    }


} // end class

